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ABSTRACT 

On the long nuclear time scale of stellar main-sequence evolution, even weak mixing processes can become relevant for redistributing 
chemical species in a star. We investigate a process of “differential heating,” which occurs when a temperature fluctuation propagates 
by radiative diffusion from the boundary of a convection zone into the adjacent radiative zone. The resulting perturbation of the 
hydrostatic equilibrium causes a flow that extends some distance from the convection zone. We study a simplified differential-heating 
problem with a static temperature fluctuation imposed on a solid boundary. The astrophysically relevant limit of a high Reynolds 
number and a low Peclet number (high thermal diffusivity) turns out to be interestingly non-intuitive. We derive a set of scaling 
relations for the stationary differential heating flow. A numerical method adapted to a high dynamic range in flow amplitude needed to 
detect weak flows is presented. Our two-dimensional simulations show that the flow reaches a stationary state and confirm the analytic 
scaling relations. These imply that the flow speed drops abruptly to a negligible value at a finite height above the source of heating. 
We approximate the mixing rate due to the differential heating flow in a star by a height-dependent diffusion coefficient and show that 
this mixing extends about 4% of the pressure scale height above the convective core of a 10 Mq zero-age main sequence star. 
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1. Introduction 

Our lack of understanding of (tnagneto)hydrodynamic transport 
processes in stars has hampered progress in developing the stel¬ 
lar evolution theory since its earliest beginnings. One particular 
aspect of the problem is the mixing in the boundary layers be¬ 
tween convection and radiative zones in stellar interiors, which is 
also known as the problem of convective overshooting. Despite 
the indisputable advance in numerical simulations, the problem 
remains extremely challenging owing to the extreme range of the 
length and time scales involved in it. 

The set of physical mechanisms that provide mixing at a con¬ 
vective/stable interface very likely depends on the type of con¬ 
vection zone involved. Because it is exposed to outer space at the 
top, a convective envelope is driven by the cold plumes originat¬ 
ing in the photosphere. It is quite possible that the plumes span 
the whole con vection zone and even pro vide mixing at its bottom 
boundary (cf. lAndrassv & SDruiil2013l and references therein). 
A convective core or shell is, on the other hand, fully embedded 
in the star, its stratification is much weaker, and the temperature 
fluctuations within it are much smaller. Consequently, a different 
set of physical mechanisms may dominate mixing at its bound¬ 
ary. 

It has long been known that the kinetic energy of the 
low-Mach-number flow in a convective core (or shell) is so 
low that the convective motions are stopped within about one 
per cent of the pressure scale height onc e they enter the 
steep entropy gradient of the radiative zone (lRoxburghlll965t 
ISaslaw & Schwarzschildlll965l) . The motions can reach much 
farther, though, if they are vigorous enough to flatten the ra¬ 
diative entropy gradient above the convective core. In this case, 
known as the process convective penetration, the motions grad¬ 
ually “erode” the radiative stratification on the thermal time 
scale until ra diative diffusion stops any further advance of the 
erosion front (IShaviv & Salneteil 1 1973b Ivan Ballegooiienlfr982t 


IZahnI 119911) . Finally, the fluid parcels hitting the stable stratifi¬ 
cation always generates a spectrum of internal gr avity waves , 
which may also provide a certain amount of mi xing (lPresslll9^ 
iGarcia Fonez & SDruidll99'TllSchatzmanlll99^ . 

Several of the processes mentioned may operate at the con¬ 
vection zone’s boundary at the same time. Their effects on long 
time scales and at long distances from the boundary are very 
different. In full numerical hydrodynamic simulations, the re¬ 
strictions on time scales that can be covered makes it difficult to 
disentangle these effects. Physical insight developed by different 
means is needed to extrapolate them to longer time scales and 
distances. 

We take a closer look at one specific process operating at a 
convective/stable interface in the interior of a star. Thermal diffu¬ 
sion causes temperature fluctuations from the convection zone’s 
boundary to spread into the stable stratification. Temperature dif¬ 
ferences on surfaces of constant pressure set up a flow even in 
the absence of momentum transport by hydrodynamic stress. We 
call this process “differential heating”, explore the physics of it 
in an idealised set-up, and estimate what amount of mixing it 
could cause in the stellar interior. 

2. The differential heating problem 

2.1. Problem formulation and simplification 

Consider a horizontal, solid surface with a stably-stratified fluid 
overlying itQ A temperature fluctuation imposed at the surface 
propagates into the fluid by a diffusive process and upsets the 
hydrostatic equilibrium. We investigate what the properties of 
the resulting flow are. 

’ Equivalently, the stably stratified fluid could be placed under the 
differentially heated surface. The role of hot and cold spots on the sur¬ 
face would be reversed in this case. We discuss only one case for the 
sake of concreteness. 
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By replacing the convective/stable interface by a solid wall, 
we eliminate all the phenomena related to the inertia of the con¬ 
vective flows and the shear induced by them. This allows us to 
study the physics of differential heating in isolation. The upper 
boundary is taken far enough not to influence the flow. Next we 
introduce further assumptions to facilitate the mathematical de¬ 
scription and the subsequent analysis of the problem: 

(1) The flow is confined to a layer that is significantly thinner 
than the pressure scale height. 

(2) The fluid is a chemically homogeneous, ideal gas. 

(3) The Brunt-Vaisala (buoyancy) frequency of the stratification 
is constant. 

(4) Thermal diffusivity is constant. 

(5) The gravitational held is homogeneous. 

(6) The differentially heated surface is flat and horizontal. 

(7) The flow is constrained to two spatial dimensions. 


Assumption (1) allows us to use the Bussinesq approxima¬ 
tion and turns out to be justified. The chemical homogeneity that 
we assume in (2) is, at least for the nuclear-burning layers in a 
star, only realistic at the onset of the burning. The differential 
heating process above a convective core would weaken as the 
nuclear burning progresses owing to the increase in the mean 
molecular weight in the core. We focus on the chemically ho¬ 
mogeneous case to keep the number of parameters tractable. We 
introduce (3) and (4) for the same reason. The Brunt-Vaisala 
frequency depends on the distance from the convective/stable 
interface in a real star. The constant frequency in our analysis 
can be thought of as a typical value for the layer influenced by 
differential heating. Finally, we add the last three assumptions 
to make our analysis more transparent and to reduce the com¬ 
putational costs of the numerical solutions. We have to keep in 
mind, however, that the constraint (7) might influence the stabil¬ 
ity properties of the flow, and thus some of our conclusions may 
not apply to the three-dimensional case. 

The Bussinesq equations are (ISniegel & Veroniill96C)t) 


V ■ u 
D16 

Dr 

'W 


0 , 

1 T' 

-Vp' H- gk + V V^it, 

Pm Tjyi 


TmN^ 


w + X V^T', 


( 1 ) 

( 2 ) 

(3) 


where u is the fluid velocity, D/Df = (9/(9fH-M V the Lagrangian 
time derivative, pm and Tu, are the mean density and tempera¬ 
ture, respectively, p' and T' the pressure and temperature per¬ 
turbations, respectively, g is the gravitational acceleration, k a 
unit vector pointed in the vertical direction, v the kinematic vis¬ 
cosity, N the Brunt-Vaisala frequency, w the vertical velocity 
component, and x the thermal diffusivity. 

Equations [T] |2] and [3] still contain several dimensional pa¬ 
rameters. It is crucial to realise that there is a natural system of 
units for the differential heating problem that makes the equa¬ 
tions dimensionless. The flow in this problem is set off by ther¬ 
mal diffusion in a stably stratified medium, therefore the inverse 
of the Brunt-Vaisala frequency, l/N (or a multiple of it), is a 
natural unit of time. Having made this choice, we can define a 
natural unit of distance as ■'/x/N, which is a typical thermal- 
diffusion length scale on the time scale 1 /N. The dimensionless 


Bussinesq equations are then 


V ■ M = 0, 

(4) 

Dm , 

- = -V» -H ?7fc + Pr V^M, 

Df ^ 

(5) 

D?7 , 

- = —W 4- V^-&, 

Df 

(6) 


where we omit any symbol to indicate the new units. We have 
also introduced a new pressure-like variable p = p'/Pm and the 
buoyancy acceleration j? = gT' jT^a, which we continue to call 
the “temperature fluctuation” in the rest of the paper, because 
that is the central concept in the differential heating process. The 
Prandtl number Pr = vfx now becomes a measure of kinematic 
viscosity, because the new unit of diffusivity is x. Equations|4l|5] 
and^are particularly well suited to theoretical studies since their 
solution is fully determined by the Prandtl number, the initial, 
and the boundary conditions. 

The distance unit ■'/xJN is rather short in stellar interiors, 
and it only weakly depends on the stratification. To see this, we 
express the Brunt-Vaisala frequency in terms of the more com¬ 
mon stellar-structure parameters, 

= ^(Vad - V), (7) 

rip 

where //p is the pressure scale height, Vad the adiabatic temper¬ 
ature gradient, and V the actual temperature gradient. Close to a 
convection zone’s boundary, we can write 


Vad - V = (8) 

where a ^ 10 ' is a coefficient of proportionality and z the dis¬ 
tance from the boundary (z > 0 in the stable stratification). When 
using Eqs.|7]and|8] the unit of distance can be expressed as 



HpHp 


1/4 


(9) 


which is about 10^ cm for values typical of a point close to the 
convective/stable interface (z IQ^^Hp) in the core of a mas¬ 
sive (IOMq), main-sequence star (x ^ 10'°cm^s“', a « 10^', 
g 10^ cm s“^, Hp ^ 10**’ cm). 

Two distinct regimes of differential heating can be expected, 
depending on the amplitude and the spatial scale of the temper¬ 
ature fluctuation imposed on the differentially heated surface. 
If the heating is strong enough, the heat transport is advection- 
dominated (i.e. the flow’s Peclet number is high), and the flow 
is generally unsteady. A similar phenomenon takes place right 
at the point where the convective flow leaves the unstable strat¬ 
ification, still retaining some positive temperature fluctuation. It 
quickly cools down as it rises in the stable medium, its tempera¬ 
ture fluctuation turns negative, and the flow is brought to a halt. 
This is the place where we can impose a lower boundary con¬ 
dition for a much weaker kind of differential-heating-induced 
flow, in which diffusive heat transport plays a major role (i.e. 
the flow’s Peclet number is low). The latter case is the main fo¬ 
cus of this paper. We show in Sect. |3] that such a flow is gen¬ 
erally smooth and reaches a stationary state (to be specified in 
Sect. 13.1b even at rather high values of the Reynolds number, 
up to Re = 4 X 10^. This allows us to gain some insight into 
the problem by exploring the scaling properties of the stationary 
differential-heating equations, which we do in the next section. 
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2.2. Analytical considerations 

The stationary differential-heating problem is described in two 
dimensions by the set of equations (cf. Eqs.|4l|5]|6]l 
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dw 


dx 

dz 
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d(uw) 
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dz 

d(uw) 


d(ww) 

dx 


dz 

d(ud) 


d{wff) 
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dp „ (d^u d^u\ 


dz 


-H T? H- Pr 


/ d^w d^w \ 

\ dx^ dz} / ’ 


d^d d^d 
dx^ dz} ’ 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 


where x and z are the horizontal and vertical coordinates, respec¬ 
tively, with the z axis pointed against the gravitational acceler¬ 
ation vector, M is the horizontal velocity component, and w the 
vertical one. In what follows, we show how the characteristic 
properties of the stationary flow depend on the typical ampli¬ 
tude 0 and the typical horizontal length scale L of the heating 
applied. 

Assume that there is a well-defined vertical length scale H 
in the differential heating flow pattern. Let us denote the typical 
horizontal and vertical velocities by U and W, respectively, and 
the typical pressure fluctuation by P. We then introduce a new 
set of variables x, z, u, w, p, and d, which all reach values of the 
order of unity close to the differentially heated surface, and 


The vertical momentum equation (Eq.[T2]i, with the substitutions 
defined above and Eqs.|2T]and|22 becomes 

d(uw) ^ d(ww) L I dp ^ ^ -I- ^ 

dx dz // \ dz / Rcjt dx^ Re^ dz} 

(24) 

and implies a close balance between the vertical component of 
the pressure gradient and the buoyancy-acceleration term pro¬ 
vided that L » // in addition to Re^ » 1 and Re, » 1. This 
allows us to estimate 

^ HQ, (25) 

which is a plain, order-of-magnitude equality of the character¬ 
istic kinetic and potential energies. Einally, the energy equation 
(Eq.fTSI) becomes 

diiid) d(wd) H „ I d^d 1 d^ 
dx dz ~~q'^^ ULdx^^WH^' 

The diffusion terms in Eq.|26]are of the order of 1 /Pex = 1 l{UL) 
and 1/Pe, = lf{WH), where Pe^^ and Pe^ are Peclet-like numbers 
associated with horizontal and vertical motions, respectively. We 
introduce them for the very same reason as we did in the case of 
Rex and Re^. Making use of Eqs.|2T]and|25] we can put Ea.l26l 
into the form 

djud) djwd) _ I d^ 1 / 

dx dz Pcx dx^ Pe^ \ L©*/^ g^} / ’ 


II 

(14) 

II 

(15) 

u — Uii, 

(16) 

w - Ww, 

(17) 

P = Pp^ 

(18) 

0 

II 

(19) 


Upon making these substitutions in Eq.[T0l we obtain 


which can be greatly simplified in the double limit of Pe^ » Pe^ 
and Pe^ « 1. In that case, the two terms in the parentheses on 
the right-hand side have to closely balance one another, so that 
we can estimate 

H ^ (28) 

and Eg.lZTlbecomes linear. 


U du W dw 
L dx H dz 

which implies the approximate relation 

— ~ — 


( 20 ) 


( 21 ) 


The horizontal momentum equation (Eq.fTTI) attains the form 


d(iiu) d(uw) P dp Pr dru Pr d^ii 

-1-~-— H-1-(22) 

dx dz dx UL dx^ WH dz} ’ 

where Eq.|2T]has been used, so the equality is only approximate. 
Nonetheless, we can see that the viscous terms are of the order 
of 1/Rex = Pr/(t/L) and l/Re^. s PiKWH), where Re^ and 
Re^ are Reynolds-like numbers associated with horizontal and 
vertical motions, respectively. We introduce this unusual nota¬ 
tion to characterise the relative contributions of the two viscous 
terms in the case of L » //. We focus on this limit because it 
turns out to be the relevant one in stellar interiors (see Sect. [S]). 
Prom now on, we assume Rex » 1 and Re^ » 1. Equation |22] 
shows that pressure fluctuations are of the order of in this 
high-Reynolds-number limit, so that we can estimate 

P ^ U^. (23) 


Equation[29lis a special case of the energy equat ion in the low- 
Peclet-number approximation of lLignie^ (Il999h . 

Using Eq.|^ we eliminate H from Eq.|25]to get an estimate 
of U{Q,L) and, with Eq.|2T] also an estimate of W{Q,L). The 
resulting relations also enable us to express Rex, Re^, Pex, and 
Pej as functions of ©, L, and Pr. This way we obtain 


U ^ 

(30) 

W Q^Pp-P} 

(31) 

Rex ^ Q'^Pp^Ppr-} 

(32) 

Re, 

(33) 

Pex Q'^Pp^P, 

(34) 

Pe^ Q^Pp-^P. 

(35) 


One might be tempted to estimate the time scale t of flow 
acceleration towards the stationary state directly from the buoy¬ 
ancy acceleration © provided by the temperature fluctuation im¬ 
posed on the bottom boundary. It is crucial to realise that, as 
Eq. |24] shows, the buoyancy acceleration is almost completely 
compensated for by the vertical component of the pressure gradi¬ 
ent in the case L » //. It is only their difference that contributes 
to the vertical acceleration. We can, however, consider the hori¬ 
zontal acceleration provided by the horizontal component of the 
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Id. 

0 

L 


Pec 

srOO 

10“ 

10* 

8.5 X 10“ 

2.5 X 10“ 

srOl 

10“ 

102 

1.4 X 102 

1.3 X 10“ 

sr02 

10“ 

102 

2.0 X 10“ 

5.5x 10-* 

sr03 

10“ 

10^ 

2.9 X 10“ 

2.7x 10-* 

srlO 

10-‘ 

10* 

2.5 X 10“ 

3.1x 10-* 

srll 

10-‘ 

102 

4.0 X 10* 

1.4x 10-* 

srl2 

10-‘ 

102 

5.5 X 102 

6.8 X 10-2 

srl3 

10-‘ 

10“ 

7.7 X 10“ 

3.5 X 10-2 

sr20 

10-2 

10* 

7.2 X 10-* 

4.1 X 10-2 

sr21 

10-2 

102 

1.1 X 10* 

1.8 X 10-2 

sr22 

10-2 

102 

1.5x 102 

9.3 X 10-2 

sr23 

10-2 

10“ 

2.1 X 10“ 

4.8 X 10-2 

sr30 

10-2 

10* 

2.0 X 10-' 

5.4 X 10-2 

sr31 

10-2 

102 

2.9 X 10“ 

2.5 X 10-2 

sr32 

10-2 

102 

4.0 X 10* 

1.3 X 10-2 

sr33 

10-2 

10“ 

5.5 X 102 

6.7 X 10-“ 


Table 1. Parameters of the series of simulations sampling a patch 
of the parameter space {0, L] at the constant value of Re = 2.6 x 

lOl 


Id. 

Resolution 

Re 

Re32 

32x32 

3.2 X 10* 

Re64 

64x 64 

6.4 X 10* 

Re 128 

128 X 128 

1.3x 102 

Re256 

256x256 

2.6 X 102 

Re512 

512x512 

5.1 X 102 

Re 1024 

1024 X 1024 

1.0 X 102 

Re2048 

2048 X 2048 

2.0 X 102 

Re4096 

4096 X 4096 

4.1 X 102 


Table 2. Parameters of the series of simulations with Re increas¬ 
ing at the fixed values of 0 = 10“^ and L — 10®. 


pressure gradient and write tZ/r a; PjL » U^jL (see Eq. |23]l. 
Using Eq.[^we obtain 

r ss 0-4/7^6/7. (36) 

Einally, we would like to point out that the characteristic 
thermal-diffusion length scale corresponding to the time scale t 
is L?!'', which scales with 0 and L in quite a differ¬ 

ent way than H does (see Ea.l28ll. This comes about because our 
estimates take the back reaction of the flow on the temperature 
distribution into account. 

2.3. Numerical solutions 

The order-of-magnitude estimates derived in the preceding sec¬ 
tion assume that the flow is stationary and that there is a well- 
defined vertical length scale in the flow pattern. We performed a 
series of time-dependent, numerical simulations of the differen¬ 
tial heating problem to confirm these assumptions and to deter¬ 
mine how the solutions depend the Reynolds number and how 
they decrease with height. 

We have developed a specialised code dedicated to the study 
of the differential-heating problem, because the problem places 
rather high demands on the numerical scheme. Eor instance, it 
has to tackle the highly diffusive nature of the flow and its high 
aspect ratio and resolve a wide dynamic range within a single 
simulation box. The code is of the finite-difference type, and 


it solves the differential-heating equations on a collocated grid 
using a variant of the MacCormack integration scheme. The 
Poisson equation for pressure, which can be derived from Eqs.|4] 
and |5] (or [37] see below), is solved by a spectral method. Heat- 
diffusion terms are treated implicitly, again by a spectral method. 
In what follows, we discuss a few selected issues related to the 
numerical solution of the differential-heating equations that need 
to be borne in mind when interpreting our results. The reader 
interested in the details of the numerical scheme is referred to 
App.ia 

We use periodic boundaries in the horizontal direction and 
force the shear stress and the vertical velocity component to 
vanish at the lower and upper boundaries of the computational 
domain. One could also use non-slip boundaries, but these are 
hardly more akin to the physical reality that motivated this study 
in the first place, so we omit this case. We impose a tempera¬ 
ture fluctuation in the form )?(x, 0) = 0 sin(7rx/L) at the bottom 
boundary and force the temperature fluctuation to vanish at the 
upper boundary. The parameters 0 and L can be identified with 
the same symbols as introduced in Sect. 12.21 

The high thermal diffusivity in the differential-heating prob¬ 
lem forces us to use long implicit time steps for the heat- 
diffusion terms, which might have an adverse effect on the ac¬ 
curacy of the results. To show that this is not the case, we re¬ 
computed the simulations sr03, sr30, and Rel024 (see Tables [1] 
and |2| and Sect. |3]l, decreasing the time step by a factor of ten. 
This brings about a change in the velocity held, which is of the 
order of 0.1% in the cases sr03 and sr30 and of the order of 
1% in the case of Rel024 (measured well away from the field’s 
zeroes). The reason for this ins ensitivity to the t ime step is the 
low Peclet number of the flow. iLigniere^ (Il999h shows that in 
the low-Pe regime, the energy equation can be approximated by 
a Poisson equation for the temperature fluctuation with w as a 
source term (see also our Eq|29ll. We do not use this approxima¬ 
tion to make our code more versatile; instead, we naturally ob¬ 
tain a close equilibrium between the terms and w in Eq.|6] 
when the Peclet number is low. This equilibrium is reached so 
quickly that details of the evolution of -d- towards the equilibrium 
become irrelevant. 

It is a well-known fact that any numerical advection scheme 
either requires adding a so-called artificial-viscosity term to 
guarantee stability or it involves some viscous behaviour implic¬ 
itly. In either case, the effective Reynold number does not even 
come close to the astrophysical regime with current computing 
facilities. The artificial viscosity (be it explicit or implicit) thus 
exceeds the physical one by a wide margin, so it demands special 
attention. 

Suppose we include an explicit viscous term as in Eq. [5] 
to model the artificial viscosity. Equations and [^ show 
that for Lx 10^ (equivalent to ~ Hp in the astrophysical 
case mentioned in Sect. l24T i we have Re, » 10 ‘^Re^ as a 
consequence of // «c L. Using equidistant grids with up to 
10^ grid points in each direction, we can achieve Re^ a; 10^. 
It follows that Re^ < 10 ' and the vertical momentum trans¬ 
port is dominated by the artificial-viscosity term. A value 
Re^ » 1 is, however, expected in stellar interiors. We use a 
simple workaround, replacing the viscous term PrV^tt by the 
anisotropic form Pr;i d^ujdx^ + d^ujdz^. The coefficients Pr;^ 
and Pr^ are re-computed at each time step from the relations 
Pix - hx max |m| Reg^y and Pr^ = li^ max Iwl Regi-y, where hx and 
are the horizontal and vertical grid spacings, respectively, and 
Regrid is the Reynolds number on the grid scale. We performed a 
few numerical tests of the code on a convection problem to deter- 
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mine that the value RegUd = 4 is a conservative trade-off between 
the amount of viscous dissipation and the code’s stability, so we 
use this value in all the simulations presented here. 

The anisotropic form of artificial viscosity enables us to 
reach Re;t » 1 and Re^ » 1 at the same time on a grid of 
reasonable size. We show in Sect.[2that the solutions with con¬ 
stant values of Pr^^ and Pr^. decay exponentially with height. The 
effective, local Reynolds number decreases in proportion to the 
flow speed, and the solutions quickly become dominated by the 
artificial viscosity. This would also happen in an (otherwise ide¬ 
alised) stellar interior at some point but that point would be much 
farther from the convection zone’s boundary. Therefore, we gen¬ 
eralise the artificial-viscosity terms, and the momentum equation 
(Eq. |5]l in 2D becomes 


Dm 

d 

du 

d 

du 


-\p + '&k — 

ox 

Prx(z)^ 

ox 

1 

|co 

+ 



where are Prx(z) and Prj(z) are proportional to e with rj being 
an adjustable parameter. We setTj = 0 when we are not interested 
in the precise vertical profiles and use 77 > 0 to suppress the 
viscous terms when examining how the solutions decrease with 
height. The latter case, 77 > 0, is a rather touchy problem because 
one has to increase 77 in a few steps, always using the (almost) 
stationary flow from a previous run as an initial condition for 
the next run. Overestimating the value of 77 can lead to a lack of 
viscous dissipation in some parts of the computational domain 
and a numerical instability ensues. Finally, the very goal that 
we want to achieve by this treatment, i.e. the flow dynamics’ 
being dominated by inertial terms at great heights, becomes an 
issue since such a flow evolves on the extremely long time scale 
corresponding to its low speed. 


3. Results 

3.1. The stationarity and structure of the flow 

Our numerical investigation of the diffusion-dominated differ¬ 
ential heating problem reveals that the flow reaches a station¬ 
ary state at all values of the Reynold number that we have been 
able to achieve (up to Re = 4 x 10^). We use the rate of change 
of the quantity Umax = max|M| (taken over the whole simula¬ 
tion box) as a convergence monitor and a basis of our criterion 
for deciding the flow’s stationarity. We show in Sect. 12.21 that 
the relevant dynamical time scale near the differentially heated 
boundary should be close to t given by Eg. (confirmed a pos¬ 
teriori, see below). We pronounce the flow stationary and stop 
the simulation once the condition 

< 10^^ (38) 

has been met at least for one time scale t. A direct implementa¬ 
tion of this condition would involve extrapolation from the time 
scale of one time step. At, to a much longer time scale t, which 
would amplify the round-off noise by a factor of r/Af » 1. 
Instead, we approximate Ea.[38]bv 

< 10^^ (39) 

where is the Euler-time-stepped solution of the equation 
dujnaxidt — (Umax “ Umax)IT- Thus, is a Smoothed version of 

Mmax, lagging behind it approximately by t in time. 

The flow in all of the simulations presented in this paper is 
composed of several layers of overturning cells with flow speed 


1 Wmax Wn 


1 uMmax 
Urytnx Ot 


rapidly decreasing from one layer to the next (see Figs. |4] or 
ISll. We characterise the flow properties close to the differentially 
heated surface by a vertical length scale H, defined as the height 
above the hottest spot at which the flow first turns over (i.e. 
h’(L/2, H) = 0), and the typical horizontal and vertical veloc¬ 
ity components U - j max(M) and W - ^ max(w), respectively, 
where the maxima are taken over the whole simulation box. The 
symbols H, U, and W can be identified with the same symbols 
as used in Sect. 12.21 The flow always reaches its maximal hor¬ 
izontal speed at the bottom boundary and the maximal vertical 
speed in the first overturning cell above the hot spot. The flow 
pattern is asymmetric, with the maximum downward flow speed 
(reached above the cold spot), max(-w), always lower than the 
maximal upward flow speed, max(w). We define the characteris¬ 
tic numbers Re^, Rej., Pe^, and Pe^ in an analogous way to what 
is used in Sect. 12.21 with the difference that now we have two 
Prandtl-like numbers Pr^ and Pr^ instead of one Prandtl number 
Pr. 

3.2. Scaling relations 

We computed a grid of 16 simulations to verify the ana¬ 
lytical relations derived in Sect. 12.21 All of these simula¬ 
tions, summarised in Table [1] have a resolution of 256 x512, 
and the vertical grid spacing was adjusted so as to obtain 
Re;t = Re^ s Re = 2.6 x 10^. The decision to fix the value of Re 
is motivated by the fact that the flow pattern turns out to be 
scalable over a large part of the parameter space provided that 
Re = const.. In other words, while changing the heating param¬ 
eters 0 and L at Re = const, does change the amplitude and 
the vertical scale of the flow, the structure of the flow, as seen 
in a system of normalised coordinates xfL and zlH, remains un¬ 
changed (see Fig. nil. We show in Fig.[T]our numerical results as 
compared with the scaling relations fitted to all but the four data 
points at L = 10*. The excluded data points do not comply well 
with the premise L'3> H and are thus expected not to follow the 
scaling relations. Allowing only the constants of proportionality 
to change in the fitting process, we obtain 


1.3 0 '^’L^^\ 

(40) 


(41) 

W = 2.7 

(42) 

t = 0.76 

(43) 


The unexpectedly good fit is a result of the flow’s scalability. 

The constants of proportionality in Eqs. I40l-l43l as well as 
the structure of the flow, depend on Re. We illustrate this de¬ 
pendence by computing a series of simulations with resolution 
ranging from 32^ to 4096^. This way, we cover about two orders 
of magnitude in the Reynolds number from 3 x 10* to4x 10^ 
(again with Re^ = Re^ s Re). We perform this experiment at 
L - const, and 0 = const., so any change in Re reflects a change 
in Pr^: and Pr^. Nevertheless, we present the dependence on Re, 
because the scalability of the flow shows that the absolute values 
of Pr;c and Pr^ do not matter. Ideally, we should choose the heat¬ 
ing parameters so as to have 0 <sc 1 and L » 1 as the scaling 
relations hold true in this limit (see Sect. 12.2b . 

Equation|43]shows, however, that the flow’s dynamical time 
scale becomes extremely long in the same limit, thus making 
any high-resolution computation unfeasible. Therefore we use 
0 = 10“^ and L - 10°, which still keeps the energy equation 
approximately linear (since Pe^^ a; Pe^ a; 10 “^), but we forgo 
having L » //. Nevertheless, we expect the changes in the flow 


5 













R. Andrassy, H. C. Spruit; Overshooting by differential heating 










Fig. 1. Dependence of the global flow characteristics on the heating amplitude 0 and length scale L. Circles show the values derived 
from numerical simulations (Table [TJ. Solid lines show the scaling relations (Eqs.|28]|30][3T] andl36ll. normalised to fit all but the 
four simulations at L = 10^ which are expected to deviate from the scaling relations. 
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Fig. 2. Dependence of the maximal horizontal velocity on the 
Reynolds number. Simulation data (circles) are connected by the 
solid line to guide the eye. The scaling law Umax is 

shown by the dashed line for comparison. 


Fig. 5. Decline with height of the relative rms vertical velocity 
in four simulations with very different heating parameters. The 
flows are reasonably close to a stationary state only up to z/H = 
2.5. 



Fig. 3. Dependence on the Reynold number of the flow’s asym¬ 
metry, characterised by the ratio of the maximum upward and 
downward flow speeds. 


with increasing Re in this case to be similar to those that would 
be seen in a simulation with L » // because of the energy equa¬ 
tion’s being linear in both cases. All of this series of simulations, 
summarised in Table |2] reach the stationary state as defined by 
Eq. Figure |2] shows that the maximum horizontal velocity 
in the computational domain slowly increases in proportion to 
j^gO.054 jjj jjjg high-Re regime. The flow also becomes increas¬ 
ingly asymmetric, as shown by the ratio of the maximum up¬ 
ward and downward flow speeds plotted as a function of Re in 
Fig. [3] The seemingly asymptotic trend changes at the highest 
Reynolds number considered, but we do not know the reason for 
this change. 

3.3. Flow at great heights 

The flow speed in all our simulations quickly decreases with 
height. Figure|5]compares the vertical profiles of the root-mean- 
square (rms; computed in the x direction) vertical velocity com¬ 


ponent, Wrms(2), in four simulations with widely disparate heat¬ 
ing parameters (srOO, sr03, sr30, and sr33). We find that Wims de¬ 
creases approximately as ^ global sense with /?„ = 1.5 

almost independently of 0 and L. We use the values H{@, L) 
given by Eq.|40]instead of those measured in the simulations to 
normalise the z coordinate, because this brings the slopes much 
closer to one another. We have to keep in mind, though, that 
these flows are reasonably close to a stationary state only up to 
zjH = 2.5, because our convergence criterion (Eq. l39l l is igno¬ 
rant of the weak flow in the upper part of the simulation box, and 
consequently, that part of the flow is still slowly evolving when 
the computation is stopped. 

The simulations discussed so far use constant artificial- 
viscosity parameters Pr^ and Pr^, which leads to a rapid decrease 
in the local Reynolds number with height (in proportion to the 
decreasing flow speed). We computed another two simulations, 
this time with 0 = 10 and L - 10'. In the first one, we set 
Ptx = const, and Pr^ = const, (the constant-Pr case hereinafter), 
just as we have done so far. In the other one, we set Pr^ oc e 
and Pr^ oc e as described in Sect. l2.31 to keep a local version of 
the Reynolds number approximately constant (the constant-Re 
case hereinafter). We increased the slope rj from 0 in a few steps 
in order to make the ratio of the rms advection terms to the rms 
viscous terms, i.e. the local Reynolds number, as independent of 
height as possible; t] - 2 turns out to be a good compromise 
in this case. There is a large-scale, residual variation by about 
a factor of four in the local Reynolds number, because the sim¬ 
ple exponential profile of the artificial viscosity is not flexible 
enough to compensate for it. Using Fig. I^we estimate that this 
variation can change the velocities by ~ 0.1 dex at most. Since 
our usual stopping condition, Eq. cannot “sense” the weak 
flow at great heights, we judge the stationarity of the flow by 
comparing the rms values of the d/dt terms to the rms values of 
all other terms that appear in Eq. and require the former to 
be significantly smaller than the latter. This way, we obtain the 
results summarised in Figs.|6]|2l and[8] The constant-Pr flow can 
be considered stationary over the whole range shown, whereas 
the constant-Re flow is only stationary for zjH < 3.8, because 
the topmost part of that flow evolves so slowly that a global os- 
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Fig. 6. Effect of two different artificial-viscosity prescriptions on the flow structure. The constant-Pr flow (Pr^, Pr^ - const.; left 
panel) is compared with the constant-Re flow {Pr^, Pr^ oc e“^"; right panel). The vertical velocity component w is in both cases plot¬ 
ted on a split-logarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a non-linear way to aid visualisation. 


cillation develops before it has reached equilibrium (see Sect. [374] 
for details). 

Figure ^illustrates that the flow is somewhat faster at z/// > 
1 in the constant-Re case, as could be expected from the mas¬ 
sive increase in the local Reynolds number by as much as two 
orders of magnitude at z + 3. Much more interesting is, how¬ 
ever, that the overturning cells in the constant-Re case become 


apparently thinner with increasing height, hence with decreasing 
local temperature fluctuation. This observation suggests that the 
scaling relations derived in Sect. 12.2! could be used locally (see 
the dependence of // on 0 in Ea.l28ll. Another piece of evidence 
for this hypothesis is shown in Fig. [ 8 ] in which we compare the 
relative rates of decrease in PrinsCz), iinns(z), and Wrms)^)- The 
envelope of i^rms(z) can be approximated well by the function 
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z/H(0, L) 


Fig. 7. Effect of two different artificial-viscosity prescriptions on 
the rms vertical velocities. Note that the constant-Pr flow can be 
considered stationary over the whole range shown whereas the 
constant-Re flow is stationary only up to z/H = 3.8. 

t?e(z) with = 1.7 for zjH < 3. We then regard 

as an estimate of the local temperature fluctuation and rewrite 
the scaling relations for the velocity components, Eqs. and 
EH to obtain their local versions, 

U,(z) ^ (44) 

We(z) ^ (45) 

where Me(z) and We(z) are expected to be good envelope models 
of Urmsiz) and Wrms(2)- In other words, we expect Ue(z) oc 
and We(z) oc _ 5^^ Indeed, these 

scalings turn out to be correct, as shown in Eig. 0 Similarly, we 
can produce a local version of Ea.l28l 

Kz) ^ (46) 



z/H(0, L) 

Fig. 8. Comparison of the rms velocities and temperature fluctu¬ 
ations in the constant-Re case with two models approximating 
their global behaviour. Solid lines show (top), Wrms (mid¬ 
dle), and i^nns (bottom). Dashed lines show the model, in which 
Me(z) ^ We(z) oc , and i?e(^) with 

Pu - Pw — ^P&, and p^ - 1.7. Dotted lines show the im¬ 
proved model given by Eqs. EH ESI and E4] with -y - 1.3. The 
coefficients of proportionality have been adjusted for each vari¬ 
able independently. 


significantly. Any detailed study of this phenomenon is certainly 
beyond the scope of this paper, but our preliminary research sug¬ 
gests that it is unlikely to be a cumulative effect induced by inter¬ 
nal gravity waves since it (1) also occurs in very small computa¬ 
tional boxes, in which all internal-wave modes are over-damped 
by radiative diffusion, and (2) the temporal spectra of the av¬ 
erage horizontal velocity are featureless at periods significantly 
shorter than that of the shear flow oscillation. 


where h(z) is a local, height-dependent estimate of a vertical 
length scale analogous to H. As a result, we expect h(z) oc 
with Ph - \p», i.e. a slow thinning of the overturning cells with 
increasing height, similar to what we observe in Eigs. |6] El and 
[8] This seemingly innocuous phenomenon has very grave conse¬ 
quences for the flow at great heights. Instead of fading out expo¬ 
nentially, it decreases even faster (see Eig. B- We expand on this 
in Sect. 14. II and derive a better model for the flow’s decline with 
height to show that the flow speed drops dramatically above a 
certain point. 

3.4. Late-time evolution of the flow 

Having continued some of our simulations for as much as 10 "^t, 
we discover an intriguing phenomenon. At first, a horizontal 
mean shear flow develops on top of the differential-heating flow. 
Its amplitude grows, and the shear flow begins to oscillate at 
some point. Einally, the oscillation saturates at an amplitude 
ranging from ~ 10“^ to ~ 10° of the differential heating flow’s 
amplitude, depending on the parameters of the simulation. The 
oscillation’s period and development time strongly decrease with 
increasing Reynolds number. They do not seem to have an up¬ 
per limit but approach IOt at Re 10°. This phenomenon most 
likely has a physical origin because decreasing the time step by 
a factor of ten does not affect the shear flow or its behaviour 


4. Interpretation of the results 


4.1. Improving the model at great heights 


After picking up the threads of Sect. I3.3I we presently find that 
the diffusion-dominated, high-Re differential-heating flow actu¬ 
ally decreases faster than exponentially with height. To see this, 
we make use of two results from Sect. I3.3I Eirst, that the scaling 
relations derived in Sect. I2.2l have their local analogues, which 
hold within the flow (compare Eqs.EHJEOl and EU with Eas.l46l 
M and EH respectively). Second, that the envelope of ??rms(2) 
can be approximated well by the function §e(z) ^ at low 

heights, where p^ is independent of © and L. This allows us to 
write 


d In )?e P§ 
dz “ ~~H' 


(47) 


The characteristic vertical scale H is linked to the heating ampli¬ 
tude 0 by Eq. ED and is thus relevant close to the differentially 
heated surface, where the typical temperature fluctuations are of 
the order of 0. A straightforward generalisation of Eq.EZlis ob¬ 
tained by replacing H by the local, height-dependent estimate 
h{z) given by Eq.EH Upon doing so, we have 


d In ??' _ 
dz' h'iz')’ 


(48) 
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where we have introduced the new variables ??'(z) = 
z' - zjH, and h'{z) - h{z)IH. By Eas. I^andl46l we have 


h'(z') ^ 


and Ea.l48]becomes 


dint?; 

dz' 




/- 1/7 


(49) 


(50) 


where p may differ slightly from because we have used an 
order-of-magnitude relation in the last step. Equation l50l shows 
that In J?'(z) decreases with a fairly constant slope over a few or¬ 
ders of magnitude, but the slope starts to change as soon as a 
wider dynamic range is considered. Since the slope is propor¬ 
tional to Ea.ISOldescribes a runaway process. Indeed, the 

solution is 


Kiz') = 


C(0) 


1/7 , 


P 


(51) 


and vanishes at a finite height of Zq = 7/y6. The constant 77 '( 0 )O^ 
must be very close to unity as )?'(0) = t?e(O)/0 ~ 1, and we can 
simplify Eq.|5T]to obtain 




(52) 


where we have returned to the non-primed variables and intro¬ 
duced a new constant y - , which is a parameter to be 

adjusted to fit the numerical data. Using the local scaling rela¬ 
tions, Eqs.|44l|45l and|46l we derive the functional dependencies 


Ue(z) oc ^1 

7/// ’ 

(53) 

VVe(z) ^1 

_Z£f 

7/// ’ 

(54) 

h{z) or 1 

y z 

~ IH' 

(55) 


The functions -d-J^z), u^iz), and We(z) are shown in Eig. [8] The 
constants of proportionality in Eqs. |52l |53] and have been 
adjusted independently, but all three functions share the value 
y - 1.3. The good fit indicates that our line of reasoning is prob¬ 
ably correct. 

Can we conclude that the flow stops at the finite height we 
have just derived? No, since the scaling relations only work in 
the high-Re regime. Provided that Re is high close to z = 0, the 
flow speed quickly decreases according to Eas. l53] and l54l until 
Re a! 1 is achieved at some height zi < 7///y. The weak flow 
above this point is supported by viscosity and gradually vanishes 
as z ^ oo. 


4.2. Allowing for a buoyancy-frequency gradient 

So far, we have assumed that the flow occurs in a particularly 
simple type of thermal stratification — one characterised by a 
typical buoyancy frequency Atyp = const. Nevertheless, we aim 
to apply our results to the immediate vicinity of a convection 
zone, i.e. to a medium, in that the buoyancy frequency rises con¬ 
tinuously from zero to a finite value. In this section, we first show 
how to estimate the value of Ntyp in such a setting and then reap¬ 
ply the techniques developed in Sect. l4.1l to demonstrate how the 
varying buoyancy frequency affects the global flow field. 

To do this, we have to recover the dependence of all the rel¬ 
evant flow properties on Ntyp by returning to a system of phys¬ 
ical units. We recall that we use 1/Myp as a unit of time and 


(ji/Ntyp)'^^ as a unit of distance, which implies that the unit of 
velocity is and the unit of acceleration (hence of ?7) is 

We use these conversion factors throughout this sec¬ 
tion without mentioning them further. The height of the bottom¬ 
most overturning cell is by Eq. @0] 


T/ph = 1.3 


\ 2/7 


bypy 


/5>l/7r2/7 
^ph ph ’ 


(56) 


where we have introduced the index ,,ph” to indicate the use of 
physical units for quantities that are dimensionless in the rest of 
our analysis. The buoyancy frequency N is now an increasing 
function of Zph and can be approximated by Eqs.|2]and[8] 


IV(Zph) = 


[Hj 



(57) 


The overall scale of the flow pattern is given by the bottommost 
overturning cell, which is thus the most important. Therefore we 
estimate Wyp = N(//ph/2), i.e. 


Ntyp - 


r ^l/2 


2 


1/2 


and combine Eqs.|56]and|58]to obtain 


//ph = 1.4 


r^2H6\ 


1/9 








2/9 


(58) 


(59) 


where we have also expanded 0ph = gAT jT^ to emphasise the 
dependence on the imposed temperature fluctuation AT jT^. We 
use the sign = in Eg. l59land also in Eqs.|60l|^ and l62] below 
to indicate that we do not expect these estimates to be off by 
more than a few tens of percent. The dependence of //ph on the 
heating amplitude and length scale is somewhat weaker in Ea.l59l 
compared with Eq. |5^ because Eq. takes into account that 
any gain in the flow’s vertical extent brings about an increase 
in the typical buoyancy frequency, which in turn makes further 
penetration harder. This effect can also be seen when we express 
the characteristic velocity components and the flow’s dynamical 
time scale in physical units. 


t/ph = 


0.8 


xg^Hl 


a 


Tm I 



Wph = 
rph ~ 



- 2/3 



(60) 

(61) 

(62) 


where the exponents have slightly changed compared with 
Eqs.|4l]|l2l andgl] 

The spatial variation of N brings on a first-order effect, too; 
that is, the stratification offers less resistance to overturning in 
the bottom part of the flow field compared with the rest of it. We 
mimic this effect by using the flow’s excellent scaling properties 
under the assumption that the flow behaves locally as if N was 
constant. Our goal is to improve upon the envelope models of 
?7rms(2), u^msiz), and >17^8(2) derived in Sect. 14. II by taking the 
dependence of N on height into account. 

Our starting point is Eq. 03 with the difference that now we 
define i?' = )/e,ph/0ph, z' = Zph///ph and IT = /iph/Z/ph- We cau¬ 
tion the reader that i/e.ph refers to a model with N = N{z) and 
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not to a direct translation of ??e that appears in Eq.|48]to physical 
units. The local vertical length scale of the flow, h{z) given by 
Eq. |4^ can be translated to physical units directly, 



(63) 

This equation, together with Eq. |56] implies 


h' ^ 

(64) 


where N' - (see Eqs.|52]and|58]l. It is evident 

that h' diverges for N' 0^, i.e. z' 0^. This effect is purely 
artificial because the divergence occurs within the bottommost 
overturning cell of the flow, and the large-scale model we are 
developing here cannot capture such local phenomena. We ig¬ 
nore the divergence for now because only /z' ' appears in Ea.l48] 
and use the same procedure as in Sect. l4.1l to derive a generalised 
version of Eq. l50l 


dlnz?; 

dz' 




(65) 


where the parameter yS has absorbed all coefficients of the order 
of unity. Its value should still be of the order of unity, but it may 
be different in this model compared with the model developed in 
Sect. 14. II By analogy to the derivation in Sect. 14. II we can write 
the solution to Eq.|^in the form 


t^e,ph(^ph) ^ 



hJ 


( 66 ) 


where we have also returned to the non-primed quantities, and 
7 - jS>?e(0) ^ parameter of the order of unity. The typical 

velocity components and the typical vertical vertical length scale 
can be estimated using the local scaling relations, Eqs. |44l |45j 
l46l and|^ We obtain 


^e.ph(^ph) ^ 
We,ph(Zph) OC 
^2ph(^ph) ^ 


Nizph) 


Myp 

A^(Zph) 


N, 


typ 


A^(Zph) 


N, 


typ 


- 2/7 


- 6/7 


- 4/7 


2 _ 7 / ^ph 
9\//ph 

I --( 

9Uph 

I _ z( 

9\//ph 


9/7 


9/7 


9/7 


(67) 

( 68 ) 

(69) 


where an explicit dependence on N appears after the transition 
to physical units. These expressions diverge for z —> 0^ where 
A —> 0^ (see Eq. l57T i. which is just another illustration of the 
envelope models’ inability to capture local phenomena (see also 
the discussion above). The bottommost part of the flow should 
in reality behave approximately as if it was in a medium with 
N = Ntyp - const., so we can cut off the problematic part of the 
N{z) profile and use, for example, the function 


N(z)^ 


Ntyp 


for 0 ^ Zph ^ 0 77ph 
for Zph ^ 2 ^ph 


(70) 


instead of N{z) in practical calculations. Doing so makes the 
right-hand sides of Eqs. |67] |68] and converge to unity as 

Zph ^ 0+. 

Just as the results of Sect. 14. li do not mean that the flow van¬ 
ishes at a finite height, neither the results of this section mean 
that. Again, the sudden drop in the typical velocities predicted 
by Egs. l67] and IMl only signifies that the flow undergoes a tran¬ 
sition to the low-Re regime at a relatively low height. Eqs. [bT] 
and 1^ cease to be usable from that point on and the weak flow 
supported by viscosity gradually vanishes as Zph ^ oo. 


5. Application to stellar conditions 

The flow in a layer of thickness /!ph(Zph) and vertical velocity 
We,ph(Zph) at distance Zph from the boundary overturns a passive 
tracer in it on a time scale = /zph/we.ph- This suggests an 
effective diffusion coefficient D^ff ^ /iphWe.ph- Eor the first layer 
above the boundary, this is 

Deff(O) = Wph//ph. (71) 

At distance Zph, Egs. 1^ and l69l give 


Deff(Zph) - Deff(O) 


Nizpb) 

- 10/7 

1 _ Z [ 

Ntyp 


9 \Hp^] 


(72) 


where Atyp is given by Eq. |58] and we have replaced A(zph) 
in Eqs. |68] and l69] by A(Zph) given by Eq. |70] as discussed in 
Sect. 14.21 The constant y is of the order of unity but cannot 
be constrained further by our present analysis. It determines 
the maximum height Zmax.ph that the mixing process can reach, 

Zmax,ph = (9/r)^^®//ph. 

Eor a specific example, consider the boundary of the core 
convection zone in a 10 M© zero age main sequence star. This en¬ 
vironment is characterised by a = d(Vad - V)/d(zph/7/p) = 0.14, 
a thermal diffusivity x = 5.9 x 10'°cm^s“', a gravitational ac¬ 
celeration g = 1.1 X 10^ cm s“^ and a pressure scale height Hp = 
2.9 X 10'°cm. A mixing-length estimate for convection in the 
core produces temperature fluctuations ATIT^ ^ 10“® on a hor¬ 
izontal length scale Lph ~ Hp. Eouationl59lthen predicts that the 
typical vertical length scale is Hph 2 x 10* cm = 7 x 10“*//p. 
The typical vertical velocity (Eg. IMTl is Wph ~ 5 x 10' cms“*. 
These numbers imply Pe;. = {Wp^Hpi^lx » 2; i.e., the bottom 
part of the flow is located right at the transition between the 
regions of advection-dominated and diffusion-dominated heat 
transport. This is not a coincidence, because we are modelling 
the region where heat leaks from the convective eddies, allow¬ 
ing them to turn over and sink back to the convection zone. Such 
a flow has to have Pe;. ss 1. Therefore, the effective diffusivity 
close to the convection zone, Deff(O) in Eq. |72] is of the same 
order as the diffusivity of heat x. Diffusivities that are several 
orders of magnitude smaller than x can be important on the long 
nuclear time scale. The maximum height reached by the differen¬ 
tial heating process on this time scale can thus be approximated 
by Zmax,ph- Assuming y = 1 we obtain Zmax.ph ~ 4 x lO^^Hp. 

Equation (172) 1 is likely to be somewhat of an overestimate 
of the actual mixing rate of the differential-heating process. The 
layers mix on the hydrodynamic time scale in their interiors, but 
as long as they are stationary, transport of the tracer between 
layers takes place by diffu sion. As in the case of semiconvective 
layering fcf. ISDruidl2013h . this reduces the effective mixing rate 
to the geometric mean of the microscopic diffusion coefficient Kt 
of the tracer and the estimate (172) 1. 

More significantly, the picture is complicated by the time de¬ 
pendence of the convective heat source. Eor the IOM 0 example, 
only the bottommost part of the flow can approach the station¬ 
ary flow speed before the heating pattern changes because the 
dynamical time scale Tph ~ 5 x 10® s (Eq. |62] with the param¬ 
eter values stated above) is of the same order as the convec¬ 
tive overturning time scale in the core. This is likely to lead to 
some form of averaging, reducing the effective amplitude of the 
source. The level of this effect can probably be investigated with 
a time-dependent simulation. 
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6. Summary 

Various observations show that there is a need for some addi¬ 
tional mixing at the interfaces between the convective and ra¬ 
diative layers of stars. Even processes that are too weak to be 
detectable in numerical hydrodynamic simulations need to be 
considered as candidate sources of this mixing, because the nu¬ 
clear time scale on the main sequence is so much longer than the 
dynamical time scale of convection, and cumulative effects are 
likely to play an important role. 

In this work, we have investigated one such weak process, 
which we call “differential heating”. The differential heating 
process occurs when radiative diffusion transports a tempera¬ 
ture fluctuation from the boundary of a convection zone into the 
neighbouring stable stratification. The resulting perturbation of 
hydrostatic equilibrium triggers a weak flow, which may provide 
mixing up to some distance from the convection zone. We inves¬ 
tigated the flow that is driven by a static temperature fluctuation 
varying sinusoidally along the solid horizontal boundary of a sta¬ 
bly stratified, thin layer of gas. This low-Peclet number problem 
(i.e. a slow flow dominated by thermal diffusion) turns out to be 
intrinsically nonlinear, in the sense that the horizontal structure 
of the flow is asymmetric. Even for symmetric boundary condi¬ 
tions, the upflow is narrower than the downflowing part for the 
flow, and the shape of the flow pattern is nearly independent of 
the amplitude of the driving temperature perturbation. 

A few additional assumptions (Sect. |2]i allow us to describe 
the problem by a set of dimensionless equations, the solution to 
which depends (apart from the boundary and initial conditions) 
only on the Prandtl number. We analysed these differential¬ 
heating equations for their scaling properties under the assump¬ 
tion that the flow is stationary (Sect. 12.2b . An astrophysically 
interesting corner of the parameter space is characterised by 
Rei » 1, Re, » 1, Pe^t » Pe^, and Pe^ 1. (The x and z 
directions have to be distinguished because such flow has a high 
aspect ratio.) In this limit we derive a set of simple relations 
(Eqs. |28] and l30l-[36ll to describe how the global flow proper¬ 
ties depend on the heating amplitude 0 and length scale L. We 
And, in particular, that the characteristic vertical length scale H 
depends only weakly on the heating parameters (Ea.l28Tl. 

We developed a dedicated numerical code to solve the equa¬ 
tions. The main difficulties are related to the highly diffusive na¬ 
ture of the flow, its high aspect ratio, and the need to resolve a 
wide dynamic range in the flow amplitude within the computa¬ 
tional box (as much as five orders of magnitude). The flow in 
our two-dimensional, time-dependent simulations reaches a sta¬ 
tionary state at all values of the Reynolds number that we have 
been able to achieve (up to Re = Re^ = Re^ = 4 x 10^). The 
flow is always composed of several layers of overturning cells, 
the shape of which depends only on the Reynolds number and 
not on the heating length scale L and amplitude 0. This property 
makes the flow scaleable in the sense that the flow field corre¬ 
sponding to some heating parameters Li, 0i can be stretched in 
space and scaled in amplitude to get a good approximation of the 
flow field corresponding to a different set of heating parameters 
L 2 , 02 provided that Re is in both cases the same. This is also 
the reason the scaling relations derived in Ea. l2.2l flt the simula¬ 
tion data remarkably well at Re = const, (see Eig.[TJ. Increasing 
the Reynolds number has little influence on the flow speed, but 
it makes the flow pattern increasingly asymmetric. 

We decrease the artificial-viscosity coefficients in the code 
with height in order to keep the Reynolds number approximately 
the same in every layer of flow cells. The numerical data show 
that the global scaling relations derived in Sect. 12.21 have their 


local analogues, which can be used within the flow. The flow 
speed’s decrease with height, being locally exponential, steepens 
with the decreasing flow amplitude according to the local scaling 
relations. Based on this we derive a model of the flow’s depen¬ 
dence on height, which closely fits the numerical data over the 
whole dynamic range that we have been able to cover (as much 
as five orders of magnitude, see Eig.|8]l. The model shows that the 
flow speed drops abruptly to a negligible value at a finite height. 
The local scaling relations also allow us to generalise our results 
to the more realistic case, in which the buoyancy frequency N 
increases with height (see Sect. 14.lb . 

We illustrated the typical scales associated with the stellar 
differential-heating process with the example of the convective 
core of a 10 M© zero-age main sequence star (see Sect. |5]l. We 
approximate the mixing due to the differential-heating flow by 
an “effective” diffusion coefficient Deff, which is of the order of 
the diffusivity of heat near the convection zone and decreases 
with height according to Eq. |72] The mixing relevant for stellar 
evolution extends about 4% of the pressure scale height above 
the convection zone. 

6 . 1. Main findings of the paper 

(1) The flow has a cellular structure and reaches a stationary 
state at all values of the Reynolds number that we have been 
able to achieve (up to Re = 4 x 10^). 

(2) Both global and local properties of the flow can be described 
by a set of simple analytical relations. 

(3) The flow speed drops abruptly to a negligible value at a finite 
height above the source of heating. 

(4) The mixing relevant for stellar evolution extends about 4% 
of the pressure scale height above the convection zone of a 
IOM 0 zero-age main sequence star. 
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Appendix A: Numerical methods 

A.1. Integration scheme 


We have adapted the standard MacCormack method to suit our specific problem. In the simplest case of a one-dimensional vector q 
of conserved quantities being advected on and equidistant grid with a spacing of Ax, MacCormack’s method can be written as 


(A.l) 

Ax 

(A.2) 

o'' + 0 ^^^ 

^n+l ^ % 

Qk - 2 ’ 

(A.3) 


where < 7 ^' is the value of q at the A:-th grid point and the n-th time step, At the time step, f(q) the flux function, and we use the 
convention that any parenthesised upper index refers to a sub-step of the method instead of a time-step index. The method is linearly 
stable provided that the CFL condition At < AxjpiA) is met, where A is the Jacobian matrix of the flux vector and p{A) is the largest 
characteristic value of A. Nonlinear stability typically requires the addition of some form of artificial viscosity. MacCormack’s 
method is second-order accurate both in space and time. 

We discretise Eqs.|4]|6l andiTTlon a collocated, two-dimensional grid of M xN cells with constant cell spacing (Ax, Az). The 
two spatial dimensions and the presence of source terms in the equations forces us to significantly extend the basic MacCormack 
scheme. We begin by advecting the vector of variables q = (m, w, 1 ?) in both spatial directions using Strang splitting. 
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where we have written out the explicit form of the flux terms. The indices k and I refer to the position along the x and z axes, 
respectively. We proceed by adding the source terms to the momentum equations. 
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(A.7) 

(A.8) 


where we use second-order-accurate central differences to keep up with the order of accuracy of the advection scheme, v/ = Prj:(z/) 
and Pi = Prj(z/) are the coefficients of our anisotropic artificial-viscosity prescription (see Sect. l2.31 l. and pi+ 1/2 - {pi + pi+\)l2. The 
new velocity field is, in general, slightly divergent. We correct for this divergence by subtracting the gradient 

of a pressure-correction field, - Af V(Ap)P\ The condition V ■ = 0 leads to a Poisson equation for the pressure 

correction. 


V ■ 

V\Apf^ = (A.9) 

Since we use central differences to compute partial derivatives, the discrete form of the Laplace operator in Eq. IA.9I should be 
derived by applying the central differences twice. That would, however, lead to a sparse operator and cause odd-even-decoupling 
problems on our collocated grid. Therefore we use the standard compact Laplacian and solve the approximate pressure-correction 
equation 


(Ap4^> 


1 ,/ 


■ 2{Apf^} + {Apf^ 


kA\,l 


(Ap)' 


( 1 ) 

k,l-l 


2(Ap)<|; + {Ap) 


( 1 ) 

k,l+\ 


(Ax)2 


(Az)2 


1 

At 


“it+i,/ “/t-i,/ 


(Id) 

4,/+i 


(Id) 

4,/-I 


2Ax 


2Az 


(A. 10) 


Equation lA. 101 is solved by a spectral solver, see Sect. IA.3] Having computed the pressure correction, we apply it to the velocity 
field. 
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The approximate nature of the pressure-correction equation (Ea. IA.lOl i causes V ■ to be small, but non-zero. Practical experience 
has shown that the residual divergence is negligibly small in the flows analysed in this paper provided that the boundary conditions 
are treated properly, see Sect. lA.2] We should also write point, but our numerical tests have shown that 

the residual divergence in the velocity field becomes much smaller if we set ^ so we use the latter form. The next step is to 

integrate the remaining two terms in the energy equation. We begin by adding the —w term. 




(A. 13) 


where its latest available value, has been used. The diffusion sub-step is given by the implicit equation 
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which is also solved by a spectral solver, see Sect. lA.3l We have thus completed the first step of the MacCormack scheme, analogous 
to Eg. lA.ll and obtained the new variables p^'^\ and The second step, which we do not do not go into detail on, differs 

from the first one at two points. Eirst, advection is done using backward-space flux differencing, as in Eq. IA.2I (compare with 
Eg. lA.ll) . Second, the pressure field is updated in this step, i.e. pfj - -i- (Ap)^^^\ The final step of the MacCormack’s scheme, 

Eg. lA.3l is used in the same form, with q = (u, w, i?). We also update the pressure field in the same way, pf^^^ = j {ph + so 
that we obtain an estimate of the pressure field for the next time step. 

Einally, there is a simple way of increasing the accuracy of the scheme at a given grid resolution, which we use. The MacCormack 
method contains a built-in asymmetry: Egs. lA.ll and lA.2l show that it always starts with forward-space flux differencing and continues 
with backward-space flux differencing. The two flux-differencing methods can be reversed, obtaining a “reverse” MacCormack 
method, without decreasing the order of accuracy of the overall scheme. We compute every time step using both the “direct” and 
the “reverse” methods and use the arithmetic average of the estimates given by the two methods. 


A.2. Boundary conditions 

The treatment of boundaries is restricted by our decision to use spectral solvers, which do not allow changing the differentiation 
operators anywhere in the computation domain. We use the ghost-cell technique for this reason. The boundary conditions we 
impose on the differential-heating flow are summarised in Sect. 12.31 The periodic boundaries in the horizontal direction are trivial 
to implement. The solid boundaries on the top and bottom of the computational domain, however, require much more care. We 
implement them using reflective boundary conditions for the velocity vector. 


Uk-\ 

= Uk,Q, 

(A. 15) 

Uk,N 

- Uk^N-l, 

(A. 16) 

Wk-l 

- -Wkfl, 

(A. 17) 

Wk,N 

= -Wk,N-l, 

(A. 18) 


so that the imaginary walls are located at I - -1/2 and at / = A - 1 /2. The conditions imposed on u also eliminate any shear on the 
boundary. The pressure field is required to be symmetric with respect to the solid boundaries. 


Pk,-i^Pk,Q, (A. 19) 

Pk,N - Pk,N-l- (A.20) 

The conditions imposed by Egs. lA.15l4A.20l can easily be shown to be consistent with the pressure-correction equation (Eg. lA.lOl 
sum both sides over k - 0, 1, ..M and / = 0, 1, ..., N). They typically do, however, bring about a cusp in the pressure field 
along the normal to the walls. The resulting discontinuity in the vertical pressure gradient then propagates to the rest of the domain 
and can be seen as a low-amplitude oscillatory field superimposed on the true pressure field (see the left panel of Eig. lATI) . We tried 
to cure this problem by changing the discretisation of the vertical-gradient operator at the walls, so that the ghost cells would not be 
used when computing the pressure gradient. This solution has met with very little success, most likely because the abrupt change 
in the operator brings about an abrupt change in the discretisation error so the problem remains. Quite surprisingly, preceding the 
pressure-gradient computation by high-order pressure extrapolation to the ghost cells has turned out to be an effective solution, 
able to eliminate nearly all of the spurious oscillations (see the middle and right panels of Eig. IA.il) . We therefore use sixth-order 
extrapolation in the simulations with constant artificial viscosity and increase the extrapolation order to ten when we let the artificial 
viscosity decrease with height. This technique cannot be viewed, however, as an all-purpose solution, because it is likely to be too 
unstable to be useful when computing highly turbulent flows. 

We require the temperature fluctuation to have a fixed sinusoidal profile at the bottom boundary and to vanish at the upper 
boundary, which translates into 
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A.3. Spectral solvers 


We use spectral methods to solve the two equations involving the Laplace operator, the Poisson equation for the pressure-correction 
equation (Eq. lA.lOb and the implicit heat-diffusion equation (Eq. IA.14b . We express both the knowns and unknowns as linear 
combinations of the Laplacian’s eigenfunctions that comply with the desired boundary conditions. The solution procedure is then 
much simplified and effective, provided that the transform to the eigenfunction basis can be computed efficiently. 

In case of the pressure-correction equation (Eg. lA. Kil l, we use the linear transform 


and its inverse 


j M-l N-1 

~ 2MN ^ ^ XI 


M-l N-l t 

" / COS 


*:=0 1=0 


7:n[l+^) 


N 


exp - 


2nimk 


M 


fkj — 'y j fm,0 "h 2 ^ j fm,n COS 
m=0 L n=l 


ni [l+{) 
N 


exp 


2nimk \ 
M ] 


(A.23) 


(A.24) 


to transform any field /jtj to an array of complex amplitudes and back. We can see that the basis functions in Eq. IA.24I are 
periodic in k and even around I - -112 and / = A - 1/2; i.e., they comply with our boundary conditions on the pressure field (see 
Sect. lA.2b . Upon using the spectral decomposition defined by Ea. lA.24l on both sides of the pressure-correction equation (Eg. I A. 1 Ob . 
we readily obtain its solution in the wavenumber space. 





(A.25) 


where we have omitted the upper indices because the expression applies to both steps of the MacCormack scheme, S k,i is the 
transformed right-hand side of Eq. lA.lOl The eigenvalues of the Laplacian are 
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(Ax)2 (Az)^ 
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and can be pre-computed. We set Aq^q to a large number to prevent division by zero and make the undetermined component (Ap)o,o 
vanish. 

In case of the heat-diffusion equation (Eq. IA.14b . we use the linear transform 
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and its inverse 
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to transform any field gkj to an array of complex amplitudes g,„^„ and back. We can see that the basis functions in Eg. lA.28l are 
periodic in k and odd around I = -1/2 and f = A - 1/2; i.e., they comply with our boundary conditions on the temperature field 
in case of a vanishing heating amplitude (see Sect. IA.2b . To allow for an arbitrary heating profile at the bottom boundary, we take 
out the known boundary term from the Laplacian on the right-hand side of Eq. IA.14I and treat it as a source term. One can show 
that it is the same as replacing the diffusion equation dd-jdt - by the equivalent equation d{§ — = V^(i7 - ^), where ( is 

the static solution to the diffusion equation d^jdt - with the desired boundary conditions (/" can be pre-computed for a fixed 
heating profile). The boundary conditions on the difference d- — ^ are then identically zero, and the spectral decomposition defined 
by Eg. IA.28l can be used. This way we obtain an explicit expression for the solution of the implicit Eg. IA.14l in the wavenumber 
space. 
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where the eigenvalues A,„ „ of the Laplacian are 


A 




2-2cos(2f^) 2-2cos(^) 

(Ajc)^ (Az)2 


(A.30) 


and can be pre-computed. An equation analogous to Eg. lA.29lrelates to 

In the practical implementation, we use the EETW library (lErigo & .lohnsonl2d()5h to compute the transforms in Egs. lA.23llA.24l 
IA.27I and lA.28l We combine standard, one-dimensional transforms of different kinds to obtain the non-standard, two-dimensional 
transforms that we need. Namely, Eq. IA.23l is implemented as a series of DCT-II transforms over the rows of the input array, after 
which the columns of the resulting array are transformed by a series of DTE transforms. The backward transform (Eg. lA.24l) is then 
computed by a series of DETs followed by a series of DCT-IIIs. The transforms for the diffusion equation (Eg. IA.27l and lA.28T l are 
implemented in the same way, but simply replacing the DCT-IIs by DST-IIs and DCT-IIIs by DST-IIIs. The transforms from the 
EETW library do not include the normalisation factor (2MA)“'. 
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Fig.A.l. Effect of three different methods of treating the pressure at the solid top and bottom boundaries. In all three panels, the 
vertical velocity component, w, is plotted on a split logarithmic colour scale. We use 0 = 10“^, L — 10\ constant kinematic viscosity 
and set the resolution to only 16 x 64 to make the spurious oscillations visible. We obtain the result plotted in the left panel using the 
simple symmetry conditions for pressure (Eqs. lA. l^ and lA.20l i. Preceding the pressure-gradient computation by third-order pressure 
extrapolation to the ghost cells reduces the oscillations’ amplitude by a factor of ~ 100 (middle panel). Increasing the extrapolation 
order to six brings about another decrease by a factor of ~ 30 in the oscillations’ amplitude (right panel). The pressure gradient is 
in all three cases computed by the second-order central differences in the whole computational domain. 
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